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ROBUST AIRFOIL OPTIMIZATION TO ACHIEVE CONSISTENT DRAG REDUCTION 

OVER A MACH RANGE * 


WU Lit, LUC HUYSEt, AND SHARON PADULA§ 

Abstract. We prove mathematically that in order to avoid point-optimization at the sampled design 
points for multipoint airfoil optimization, the number of design points must be greater than the number of 
free-design variables. To overcome point-optimization at the sampled design points, a robust airfoil optimiza- 
tion method (called the profile optimization method) is developed and analyzed. This optimization method 
aims at a consistent drag reduction over a given Mach range and has three advantages: (a) it prevents severe 
degradation in the off-design performance by using a smart descent direction in each optimization iteration, 
(b) there is no random airfoil shape distortion for any iterate it generates, and (c) it allows a designer to make 
a trade-off between a truly optimized airfoil and the amount of computing time consumed. For illustration 
purposes, we use the profile optimization method to solve a lift-constrained drag minimization problem for 
2-D airfoil in Euler flow with 20 free-design variables. A comparison with other airfoil optimization methods 
is also included. 

Key words, robust optimization, airfoil shape optimization, consistent drag reduction, lift-constrained 
drag minimization, adaptive weight adjustment 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. Optimization of aerospace vehicles or aircraft wings is based on a mathematical 
model of the physical reality. The design variables are parameters in the mathematical model and changes 
in these design variables result in new physical objects. Aerodynamic optimization is a process of finding 
a set of design variables that corresponds to a new physical object with better aerodynamic and structural 
properties. 

For airfoil shape optimization, design variables are parameters that define the airfoil shape. One accepted 
practice for modeling airfoil shapes is to use empirical algebraic expressions that are based on the knowledge 
of aerodynamic properties of airfoils. There are two advantages to such an airfoil model: 

(a) Each set of design variables (such as maximum thickness, camber, radius of leading edge, etc.) 
generates an airfoil with aerodynamic properties that are well understood by experts. 

(b) The number of design variables is small, thus the corresponding optimization problem is computa- 
tionally affordable. 

A drawback to such an airfoil model is that the optimization process is merely selecting a desirable airfoil 
shape among the predetermined shapes given by the specific airfoil model. The usefulness of the optimal 
airfoil is essentially determined by the usefulness of the airfoil model. 

To achieve truly innovative airfoil designs, it is therefore desirable to consider “free-form” shape opti- 
mization, which allows “truly optimum” designs to be computed [5, Section 1]. Here “free-form” means 
geometric shapes represented by a linear combination of general basis functions (such as splines or sinusoidal 
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basis functions). However, there are several challenges associated with “free-form” airfoil shape optimiza- 
tion. Drela [5, Subsection 5.1] pointed out that if presented with sufficient design model resolution, an 
optimizer will readily (and annoyingly) manipulate and exploit the flow at the smallest significant physical 
scales present that tends to produce improved performance only near the sampled operating conditions. 
The point-optimized airfoil often shows a possibly severe degradation in the off-design performance and 
optimized aerodynamic shapes are usually “noisy” and usually require a posteriori smoothing. 

The main objective of this paper is to develop a robust airfoil optimization scheme for achieving consistent 
drag reduction over a given Mach range. This scheme (called the profile optimization method) has the 
following three advantages: (a) it prevents severe degradation in the off-design performance by using a smart 
descent direction in each optimization iteration, (b) there is no random airfoil shape distortion for any iterate 
it generates, and (c) it allows a designer to make a trade-off between a truly optimized airfoil and the amount 
of computing time consumed. 

The term “robustness” has been used to mean a variety of things. For optimization under uncertainty, 
we can distinguish the following meanings and goals of “robust optimization” : 

1. Identify designs that minimize the variability of the performance under uncertain operating condi- 
tions. This is the objective of Taguchi methods [8], which are most practical when the expected 
value of the performance can be adjusted at negligible cost. 

2. Mitigate the detrimental effects of the worst-case performance. This is the objective of minimax 
strategies, which can be used to find a design with the optimal worst-case performance [4], 

3. Provide the best overall performance of a system by maximizing the expected value of its utility 

[ 12 ]. 

4. Achieve consistent improvements of the performance over a given range of uncertainty parameters. 
This is the main objective of this work. 

The paper is organized as follows. In Section 2, we give a brief introduction to the airfoil optimization 
problem. Section 3 is devoted to Drela’s hypothesis on the necessity of using at least 0(m) design points 
for multipoint airfoil shape optimization, where m is the number of free-design variables. The main result 
is a mathematical proof of the fact that the multipoint airfoil shape optimization needs at least (m + 1) 
design points. We present two robust airfoil shape optimization formulations in Section 4 and prove the 
mathematical equivalence between these two formulations under finite sampling in Section 5. The profile 
optimization method is introduced in Section 6. The implementation issues of the profile optimization 
method are discussed in Section 7 and numerical simulation results are given in Section 8. Final conclusions 
are drawn in Section 9. 

2. Airfoil Shape Optimization. Recently, there has been significant progress in airfoil shape opti- 
mization (see [2, 3, 5, 16] and references therein). These papers demonstrate impressive shape optimization 
using high-fidelity CFD codes, reliable grid generation, and numerically efficient sensitivity calculations. 
Equally impressive progress has been made in optimization of 3-D wings [6, 7, 17] and in coupled structural- 
aerodynamic optimization [9]. However, except Drela’s work [5], these aerodynamic shape optimization 
projects all find optimal shapes based on fixed operating conditions. In this paper, we study a simplified 
2-D airfoil shape optimization problem using a low-fidelity Euler flow solver, but we include uncertainty in 
the operating conditions for the airfoil shape optimization. 

Airfoil shape optimization is a PDE-const rained optimization problem. A general mathematical frame- 
work for airfoil optimization can be described as follows (see [3, 16, 17] for details). 

For a given flow and turbulence model (a set of PDEs), one could compute the energy, velocity, and 
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pressure around a given airfoil, which will be denoted by one vector-valued function S. For an airfoil 
optimization problem (which is a 2-D problem), S is a function of (x,y) for fixed ( D,a,M ), where D 
denotes the set of geometric design parameters that defines the airfoil, a is the angle of attack, and M is 
a given Mach number. Note that M and a are parameters representing the flight conditions. The PDE 
equations for S are solved by a numerical method that involves grid generation, discretization of the PDEs, 
and iterations for solving systems of nonlinear equations. 

For a given geometric design vector D , the generated grid (denoted by X ) is determined by D, i.e., 
X — X(D) is a function of D. Usually, X(D°) is generated for an initial design vector D° based on a grid 
generation method, and X(D) (for D / D°) is generated by a grid movement strategy (such as the spring 
analogy or the elasticity model) [16]. 

The system of nonlinear equations derived from a discretization of the PDEs for the flow and turbulence 
model will be denoted by 

(2.1) R(D, a, M, X, S) — 0, 

where S is a vector whose components are values of S at the grid points (or the cell centers). For any given 
(D, a, M, X), a numerical solution of (2.1) gives S. 

If we consider minimization of the drag coefficient under a minimal requirement on the lift coefficient, 
then the multipoint airfoil optimization problem can be formulated as follows [5]: 

r 

(2.2) min ^ wi c d ( D , a u Mi,X(D), Si(D , a*)) 

i—1 

subject to 

(2.3) ci(D,ati,Mi,X(D),Si(D,ai )) > c,* for 1 < i < r 

and D € T, where iDj’s are positive weights, cf is the minimal lift value, T is a given feasible set for 
geometric design variables, Si(D,oti) is the solution of (2.1) for M — M t and a — a,, and Cd (or cf) denotes 
an approximate value of the drag (or lift) coefficient. The feasible set T mainly depends on the airfoil model. 
In this paper, we use splines for airfoil shape modeling: the components of D are the control points for a 
spline curve and the feasible set T — JR”' , where m is the number of geometric design variables. 

Recently, Drela [5] studied the behavior of the optimization solutions of (2.2) in two-dimensional viscous 
flow when the number of free-design variables is relatively large. Drela [5] concluded that increasing the 
number of geometric design variables requires a corresponding increase in the number of design conditions 
(Mach numbers) used in the multipoint optimization problem (2.2). He also suggested that it is necessary to 
have r — 0{m), where m is the number of free-design variables, to achieve a smooth airfoil geometry. Other 
notable conclusions made by Drela [5, Conclusion] are: 

• Near-continuous sampling of the operating space (i.e., in the range of Mach numbers) may be 
required in the theoretical limit of a general airfoil design problem with a very large number of 
degree of freedom (for geometric variables) — a very expensive proposition. 

• The most suitable operating points to be actually sampled in multipoint airfoil optimization (i.e., 
Mi, M 2 , . . . , M r ) are not apparent a priori. From limited experience, sampling somewhat beyond 
the expected operating range appears to be best. 

• The point weights (i.e., W\, w 2 , ■ ■ ■ , w r ) used in multipoint airfoil optimization are arbitrary, and 
their appropriate values can not be easily estimated without prior experience. 

• Optimized aerodynamic shapes are usually “noisy” and require a posteriori smoothing. 


3 



3. The Critical Number of Design Points for Multipoint Airfoil Optimization. In this sec- 
tion, we give a mathematical proof of Drela’s hypothesis on the necessity of having the number of design 
points proportional to the number of design variables. Specifically, we prove that in order to avoid point- 
optimization at the design points, the number of design points must be at least (m + 1), where m is the 
number of free-design variables. 

In (2.2), the angle of attack, a*, is predominantly determined by the constraint on the lift coefficient at 
M t . Therefore, it is theoretically possible to eliminate the constraint and consider the following unconstrained 
reformulation of (2.2): 

r 

(3.1) min Z 

i=i 

where Mj’s are design points, Wi s are positive weights, and / is a function related to the drag coefficient. 
Let D be an optimal solution of (3.1). Then 

(3-2) 0, 


where denotes the gradient of / with respect to D. 

If r < m, then (r — 1) < in. So we can find a nonzero vector A D in the orthogonal complement of the 
subspace of R m generated by the following (r — 1) vectors: 




By (3.2) and w\ > 0, AD must also be orthogonal to |jy(-D,Mi). Therefore, we obtain a nonzero vector 
AD such that 


(3.3) 


(2L(D,Mt),AD)=0 


for 1 < i < r, 


where (u, v) denotes the dot product of vectors u and v. 
Let M r+ 1 be any Mach number such that 


(3.4) 



(D,M r+1 ),AD 


^0. 


Let {w i, . . . , 'i£’ r +i } be a set of arbitrary positive weights. By Taylor expansion, we get 


r+1 


r+1 


r+1 


(3.5) 


i— 1 1=1 

It follows from (3.5) and (3.3) that 


i= 1 


Y, + tAD ’ M i ) = Z mo + M «)> AD ) + °A)- 


dD' 


r+1 


r+1 


(3.6) 


Z + tAD , Mi) = z Mi) + ^r +1 ( |£(A M r+1 ), AD ) + 0 (t 2 ). 


8D' 


By (3.4), we can choose t such that |f| is sufficiently small and 

t|| ( Atf ffl ) 1 Au) <0. 
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Then (3.6) implies 


r+1 r+1 

£ mf(D + tAD, Mi) < mf(D, Mi). 

i=l i = 1 

That is, after adding one more design point, D is not an optimal solution no matter how the weights are 
selected. ■ 

In summary, if the number of design points is less than the number of free-design variables, then it is 
possible to have a first-order reduction of the value of / at some Mach number M r+1 (in order of |t|) and only 
second-order increments of the values of / at M t ! s (1 < i <r) (in order of t 2 ). This results in a much better 
performance of the airfoil at an off-design point M r+1 with a marginal deterioration of the performance at 
other Mach numbers. 

4. Robust Airfoil Optimization Formulations. For engineering design problems, an optimal design 
is usually obtained under some explicit/implicit assumptions. This leads to a design that works well under 
ideal operating conditions but may perform inadequately under non-ideal (i.e., off-design) conditions. The 
problem is that the optimal design does not consider the uncertainty or variability of some parameters/data 
that will affect the actual performance of the design in a real-world situation. Therefore, it is necessary to 
include uncertainties in a practical design optimization process. In this paper, we assume that M is the only 
uncertain parameter and [M m j n , M max ] is the range of Mach numbers considered. 

Let p(M) be the probability density function of the uncertain parameter M. Then a stochastic pro- 
gramming formulation for airfoil optimization under uncertainty can be described as follows [12]: 

pMm ax 

(4.1) min / c d (D,M,a(M))p(M) dM 

D M) J Mmi„ 

subject to 

(4.2) c t {D, M,a(M)) > c ( * for M min < M < M msx 

and D e T, where T is the feasible geometric design space and c\ is the minimal lift value. This corresponds 
to the third definition of “robust optimization” given in the introduction. 

On the other hand, one can also use the following minimax optimization formulation for robust opti- 
mization under uncertainty [4]: 

(4.3) min max p(M) c d (D , M , a(M)) 

^min ^ ^ ^ ^max 

subject to the constraints given in (4.2). Here p(M) > 0 is a positive weighting function of M. This 
corresponds to the second definition of “robust optimization” given in the introduction. 

The constraints on the lift can be eliminated by choosing an appropriate value for the angle of attack 
corresponding to each M. In fact, Elliott and Peraire [6] suggested to incorporate a means for adjusting 
the angle of attack to satisfy the lift constraint into the flow analysis algorithm. One problem with this 
elimination approach is that the angle of attack becomes a function of M and D. and the derivatives of a 
with respect to D have to be computed in order to get the derivatives of c; and c d with respect to D. In 
this study, we retain the lift constraint to make our methods applicable to general constrained aerodynamic 
optimization problems (with constraints on the pitching moment and the leading edge radius, etc.). 
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5. Approximations to Robust Optimization Formulations. Since we cannot compute C; and c,j 
for all M in the Mach range [M m j n , M max ], computationally tractable approximations of (4.1) and (4.3) must 
be used. The simplest approximation scheme is to replace [M m j„, M max ] by a finite subset of [Af min . M max ], 
say {Mi, M 2 , ■■■ , M r }. Then (4.1) is reduced to the following multipoint optimization problem [5]: 

r 

(5.1) min Y' WiCd{D,Mi,ai) 

D,a i,...,a r 

i—1 

subject to 


(5.2) D e T, ci{D,Mi,a.i) > til for 1 < i < r, 

where the weights Wi depend on the probability density function p(M) and a chosen integration scheme. 
Similarly, the minimax formulation (4.3) can be discretized as follows: 


(5.3) min max piC d (D,Mi,ai), 

D,a i,...,av l<i<r 

subject to the constraints given in (5.2). Here pi > 0 is determined by p{M) and M , . 

It seems that (5.1) and (5.3) are completely different. However, under the strict complementarity 
condition ( i.e ., Lagrange multipliers are nonzero for all active constraints), (5.1) is mathematically equivalent 
to (5.3). Here we give a proof of the mathematical equivalence between (5.1) and (5.3). 

Let ( D , qu,..., Q r ) be a stationary solution of (5.1). For simplicity, assume that D is in the interior of T 
and the equality holds for lift constraints. Then, under the strict complementarity condition, the following 
first-order optimality condition (or the KKT-condition) for (5.1) holds: 


(5.4) 


^m^&Mu&i) -jrx^&Mu&i) - 0, 

i = 1 i—1 

~ dci 

Wi-~-{D, Mi, a,) - A i~{D, M^&i) for 1 < i < r, 

OOti C/Cki 

A i > 0, Ci{D, Mi, a i) — Ci for 1 < i < r. 


Define 

r 

7 == y^WjC d (D,M i ,dj), 

i= 1 


pi ^ >0 for 1 < i < r, 

Cd{D, Mi, &i) 

6i — for 1 < i < r. 

Pi 

Then (5.4) implies the following conditions: 

r ft,,, r „ 

X] ^pi~d£( D , «») - E x >db {D ! Mu &i) = 0> 


QiPiw^~{D,Mi,ai) — A i-—l{D,Mi,dii) for 1 < i < r, 

OOti OOti 

(5.5) A i > 0, ci(D, Mi, oti ) — Cj* for 1 < i < r, 
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piC d (D, Mi,&i ) =7 for 1 < i < r, 


8i > 0 for 1 < i < r, and ^ 9i — 1 . 

i=l 

By (5.5), (7, D, di, . . . , a r ) satisfies the KKT-condition of the following optimization problem: 

'V 

(5.6) min 7 subject to ci(D,ai,Mi) > d{ and c d (D,ai, M t ) < -L for 1 < i < r. 

D,a i,7 pi 

Since (5.6) is mathematically equivalent to (5.3), ( D , au, . . . , a r ) is also a stationary point of (5.3). 

On the other hand, assume that (£>, dq , . . . , a r ) is a stationary point of (5.3), the equality holds in (5.2), 
and piC d (D, Mi, &i) = 7 for 1 < i < r, where 

7 : = max c d (D,Mi,&i). 
l<i<r 

Then (7, D, di, . . . , d r ) is a stationary solution of (5.6). If we further assume that the strict complementarity 
condition holds, then the KKT-condition (5.5) for (5.6) holds. It follows from (5.5) that (5.4) holds for 
Wi = piOi. Therefore, (D,a.\ , ■ ■ . , a r ) is also a stationary point of (5.1). ■ 

Remark 1. Note that the strict complementarity condition holds naturally for (5.1). In fact, if the 
lift is greater than its target value, then one can reduce the angle of attack so that the lift with respect to 
the new angle of attack is equal to the target value, but the drag with respect to the new angle of attack 
is reduced. This is not possible at a stationary point. Therefore, all lift constraints at a stationary point 
become equality constraints. Also, the derivative of the drag with respect to the angle of attack is positive for 
the range of the angle of attack considered herein, which implies that the Lagrange multipliers X, are positive 
(cf. (5.4) or (5.5)). Thus, the strict complementarity condition for (5.1) always holds. An implication is 
that it is possible to recover an optimal solution of (5.1) by solving (5.3) with appropriate choices of pi’s. 

However, it is a little more difficult to claim that a stationary point of (5.3) is also a stationary point 
of (5.1). In the proof given above, we used two additional assumptions: (a) all piC d (D, M t ,a ») have the 
same value 7, and (b) the strict complementarity condition holds for the drag constraints in (5.6). The first 
assumption is not realistic in the sense that we do not know how to choose pi’s to make piC d (D, M t , 6q) have 
the same value for an optimal solution. Adaptive adjustments of pi’s are used in the profile optimization 
method to ensure all piC d (D, Mi, on) (1 <i <r) have the same value during the optimization iterations. The 
assumption (b) is a natural one, since the strict complementarity condition for the drag constraints means 
that every drag constraint is necessary for getting the optimal solution. ■ 

The results in [5] indicate that the multipoint airfoil optimization formulation (5.1) tends to produce 
airfoils that have possibly severe degradation in the off-design performance, which Drela referred to as point- 
optimization behavior. The above equivalence analysis also implies that the same conclusion holds for the 
minimax formulation (5.3). 

Huyse and Lewis [12] concluded that the point-optimization behavior can be attributed to the discretiza- 
tion error between (4.1) and (5.1). The multipoint formulation (5.1) approximates the integral as a discrete 
sum: low drag is obtained at the selected Mach numbers M\, . . . ,M r , but there is no control over, nor 
requirements for, the drag at the other Mach numbers. Due to the highly nonlinear nature of the PDEs 
for flow simulation, the optimizer is able to mold the objective function c d to its own advantage: distinct 
drag troughs appear at each of the Mach numbers M\, . . . ,M r (see also [5]). To prevent the optimization 
algorithm from exploiting the approximation error, Huyse and Lewis [12] suggest changing the integration 
points Mi , . . . , M r during each optimization step. This effectively ensures that low drag is obtained not 
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for just r fixed Mach numbers, but for any combination of r Mach numbers Mi, ... , M r . In their study, 
they added a random perturbation to the integration points in each optimization step but state that “any 
adaptive scheme that varies the location of the integration points M k for each optimization step will do.” 

We will show that with only far fewer design points than the number of free-design variables, the profile 
optimization method can adaptively adjust pi s in (5.3) to achieve a consistent drag reduction over the given 
Mach range, so there will be no degradation in the off-design performance of an optimal airfoil. 

6. Profile Optimization Method. In this section, we first give a motivation for a new shape opti- 
mization strategy based on the robust optimization model (5.3) and then present the profile optimization 
method. 

Let D be a fixed feasible design vector. We can plot the drag with respect to the Mach numbers over 
the interval [Af m j n , M max ] while keeping the lift at a constant value by adjusting the angles of attack for 
different Mach numbers. Such a plot will be called a profile for D (see Fig. 8.3(a) for example). Ideally, we 
want to obtain a design vector D with a desirable profile (for example, with the drag rise occurring at a very 
high Mach number). 

Suppose that a “perfect” airfoil D could be obtained by some magic process, which might be either 
theoretical or experimental. Let us look at the unknown process from a reverse engineering point of view. 
That is, even though we have no idea about how the airfoil was designed, we want to construct an optimization 
scheme that produces an airfoil with a similar or better profile. 

Assume that the profile curve is generated by computing the drag coefficients at Mach numbers Mi, M 2 , 
. . . , M r . Let 1/ pi be the value of the drag coefficient of the “perfect” airfoil at Mach number M* with angle 
of attack on, which is chosen so that the lift coefficient is equal to c(. Then the global optimal value of (5.3) 
is no more than 1. If (D, ai , . . . , a r ) is a global optimal solution of (5.3) , then 

c d (b,Mi,&i ) < — - c d (D,Mi,&i) for 1 < i < r. 

Pi 

Therefore, with the given analysis (i.e., using the profile of the drag coefficient at Mi , ... , M r ), D is no worse 
than the “perfect” design vector D. 

This shows that for any given profile analysis method, it is possible to recover or improve a desirable 
airfoil by solving (5.3), if the weights are appropriately chosen. In this sense, (5.3) provides a very flexible 
tool for airfoil optimization. By experimenting with various choices of the weights pi, a robust airfoil may 
be obtained (if such an airfoil exists). However, instead of guessing the appropriate weights in (5.3), we 
can adaptively adjust the weights during the optimization iterations. This leads to the profile optimization 
method. 

Before presenting the profile optimization method, we give a linear programming formulation of a trust- 
region method for solving (5.3). Using linearization of the nonlinear functions in the mathematically equiv- 
alent formulation (5.6) at the current iteration point ( D k , a\ j ; , . . . , a r , k ), we get the following trust-region 
subproblem for (5.6): 


(6.1) 


min 7 subject to 

AD,4«i,7 


— <7j S < A Di < OiS for 1 < i < m, 

-a itk < Acq < a max - a itk for 1 < i < r, 


c l {D k ,a uk ,M i ) + (^{D k , 


^i) > AD 




Oii, k ,Mi)A.ai > Cf for 1 < i < r, 
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Cd(D k ,a^k, 


f r\ \ 

Mi) + (^ j (D k ,a i , k ,M i ),AD\ + -^{D\a iM Mi)Aa 


t < — for 1 < i < r, 
Pi 


where AD and Aa* are the increments for D and a*, respectively, a t and 8 are positive constants that define 
the trust region for AD. and a max is the maximum angle of attack allowed. 

Algorithm 6.1. (Profile Optimization Method) Suppose that D° is a given design vector and 
Mi , M 2 , . . . , M r are the design points. Let 0 < rj < 1 be the predicted percentage reduction rate of the drag for 
the trust-region method and let 0 < e < 1 he a parameter for termination of the algorithm. Then construct 
a new design vector as follows. 

(6.1.0) Initialize the angles of attack: Find 01 , 0 , 0 : 2 , 0 , •••, a r ,o such that c;(Z) 0 , Oj,o, M;) = cf for 1 < 
i < r; and let k — 0. 

(6.1.1) Adjust weights: Let 


1 

Pl c d (D k ,a itk ,Mi) 


for 1 < i < r. 


(6.1.2) Check early termination: If the zero vector is an optimal solution of (6.1), then output D k as 
an optimal solution and terminate the algorithm. 

(6.1.3) Find a trust region for the predicted percentage reduction of the drag: Find 8 > 0 such 
that the optimal objective function value of (6.1) is (1 — rj). 

(6.1.4) Solve the trust-region subproblem: Find the least norm solution (AD k , Aoi,^, . . . , Aa r<k ) of 

( 6 . 1 ). 

(6.1.5) Generate the new iterate: Let Oj,j, + 1 = Oj,t + Ao*,*, for 1 < i <r, and D k+1 — D k + A D k . 

(6.1.6) Check heuristic termination conditions: If 


(6.2) max piC d (D k+1 ,cti, k+ i,Mi) > 1 and 

l<i<r 

r 

(6.3) e k < c d (D k ,<*»,*, M*), 

i=i 

where e k is the accumulative reduction of the drag at the design points defined by 

r 

:= ^ ( 'c d (D k , Oj, fc , Mf) - c d (D k+1 , Oj,fc + i , M») j , 

i=l 

then output D k as an optimal solution and terminate the algorithm. 

(6.1.7) Start a new iteration: Update k := k + 1 and go back to (6.1.1). 

Remark 2. The adaptive adjustment of weights makes it possible to achieve a consistent drag reduction 
over the Mach range [M m j n , M max ] . Note that moving along the descent direction A D k gives a simultaneous 
drag reduction at all design points Mi, . . . , M r (at least if a small step is used). If the selected design points 
are a “fair representation” of all Mach numbers in [ M m \ n , M max ], then moving along A D k will not induce any 
hidden rise of the drag at some unselected Mach number. This leads to a robust optimization method that 
achieves a consistent drag reduction over the given Mach range from iteration to iteration (see the profiles 
of the drag coefficient for iterates generated by the profile optimization method in Fig. 8.4(a)). 

In contrast, without adaptive weight adjustments, it is necessary to use at least (m + 1) design points, 
where m is the number of free-design variables, in order to avoid point-optimization at selected design points 
(cf. Sections 3 and 5). 

Besides the adaptive adjustment of weights, there are three features of the profile optimization method 
that are not standard for a trust-region method. 
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(a) First, the size of the trust region is modified to achieve a given predicted percentage reduction rate >/ 
of the drag. This approach is employed to avoid the following two problems, which are incurred in practical 
applications of a standard trust-region method: (i) an appropriate size of the trust region is unknown without 
prior experiences, and (ii) if each iteration takes a long time to complete, then rejecting a new iterate might 
be hard to accept by a designer. 

For airfoil shape optimization, if the percentage reduction of Cd is too small, then the reduction could be 
a consequence of numerical errors and is less reliable. If the predicted percentage reduction of Cd is too large, 
then the size of the trust region is also relatively large and the predicted reduction can not be trusted. By 
choosing an appropriate percentage reduction rate for each iteration, these two problems might be avoided. 

Note that the choice of the predicted percentage reduction rate depends on a user’s knowledge of the 
accuracy of the simulation analysis tool involved, the time allowed to generate a new design, and the expected 
overall performance improvement. For example, if the relative numerical error for computing Cd is 0.5%, 
each iteration takes 3 hours, a designer has 1 day to generate a new design, and the expected improvement 
is 8%, then the designer could set r) — 1% so that each iteration will produce some meaningful improvement 
of the design and the final design might have about 8% improvement over the original design after 1 day (or 
8 iterations). 

(b) The second non-standard feature is that we use the least norm solution of the linear programming 
problem (6.1), which can be obtained by solving a quadratic perturbation of (6.1) [14]. Our implicit assump- 
tion here is that the original airfoil is reasonable. Therefore, we do not want a new airfoil to deviate too much 
from the original one if unnecessary. By using the least norm solution of (6.1), we intend to select a new 
airfoil closest to the original one while achieving a predetermined amount of improvement in performance. 

(c) The third non-standard feature of the profile optimization method is the termination criterion, which 

is based on design heuristics. The basic idea is that if a good descent direction for drag reduction at all 
design points does not work (i.e., (6.2) holds) and the overall percentage reduction is insignificant ( i.e ., (6.3) 
holds), then there is no need to continue the optimization process. ■ 

7. Implementation of Profile Optimization Method. The main implementation issues related to 
the profile optimization method are the following: 

• parallel processing of the function and gradient evaluations for r Mach numbers, 

• finding initial values of the angles of attack for r Mach numbers, 

• finding the size r) for the trust region in an iteration, 

• computing the least norm solution of (6.1). 

In the next four subsections, we discuss these four implementation issues. 

7.1. Parallel Processing of Function and Gradient Evaluations. Approximate values of ci (or 
Cd) can be obtained with one flow analysis (i.e., solving the PDE once) and the gradient of ci (or c,i) can 
be obtained with one sensitivity analysis (i.e., solving one adjoint equation). The function values and the 
gradients of Cd (or c;) at different Mach numbers can be calculated independently. Therefore, if 2 r processors 
are available in a parallel computing environment, the walltime for function and gradient evaluations in each 
iteration of the profile optimization is one flow analysis and one sensitivity analysis. 

It is desirable to have a parallel wrapper code for the profile optimization method to run several instances 
of a CFD code without requiring any modifications to the code. For this purpose, we store multiple copies 
of an executable CFD code (FUN2D [15] for the current implementation) under different directories (called 
hosting directories) in a file server. Each copy of the CFD code communicates with the profile optimization 
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method by writing its output into its own hosting directory and the profile optimization method collects 
all function values and gradients from the hosting directories. The architecture of a parallel computing 
environment (such as shared memory or distributed memory) has no impact on the implementation since 
there is no communication between any two different processors. 

The current parallel wrapper code is based on a shell-level parallel computing protocol “pbsdsh” on Coral 
[13] at ICASE. Coral is a Beowulf-class cluster of 96 Intel Pentium processors with 3 additional file servers. 
To get the function values and the gradients for c; and Cd at r different Mach numbers simultaneously, we 
just create 2 r hosting directories. For 1 < it < r. in the i-th directory, we use FUN2D to compute the 
function value and the gradient of c; for the i-th Mach number Mj, and in the (r + i)-th directory, we use 
FUN2D to compute the function value and the gradient of Cd for the i-tli Mach number Mj. Therefore, with 
2 r processors, we get all function values and gradients for ci and c<j in the walltime for one flow analysis and 
one sensitivity analysis. 

Note that on Coral, the Unix command “pbsdsh gradfun” makes the same executable code “gradfun” 
run on different nodes of Coral in parallel. Here is a brief description of how to write a C-code that generates 
an executable code “gradfun” for parallel processing of function and gradient evaluations on Coral: 

• Create a shell script called “get jd” , whose content is 

echo $PBD_NODENUM 

The output of “get Jd” is the value of the environment variable $PBD_NODENUM, which is (i — 1) 
if “get Jd” is executed on the i-th node. 

• Write a C-code that does the following: 

- use “gethostname(HOST-NAME)” to get the host name (such as n013) of the node running 
the C-code; 

- execute a system command 

“getJd > HOST_NAME” 

in the C-code that gets the SPBDJMODENUM environment variable and stores it in a file 
whose name is the content of the string HOSTJNAME; 

- read the index number id in the file whose name is the content of HOST_NAME; 

- go to the id-th hosting directory; and 

- execute FUN 2D to get the function value and gradient for c; (or a) at Mj if and only if i = id+1 
(or r + i = id + 1) for 1 < i < r. 


7.2. Finding Initial Values of the Angles of Attack. For a fixed airfoil shape and a fixed Mach 
number, the lift is almost a linear function of the angle of attack [10, p. 99]. Therefore, we use the following 
search method based on a linear interpolation to find an angle of attack for which the corresponding lift is 
the given target value. 

Algorithm 7.1. Choose an error tolerance e > 0 and an initial adjustment rate 0 < k < 1. Let ao be 
an initial guess of the angle of attack for a given Mach number M and a given geometric design vector D. 
Then the following procedure will find a value of a such that \ci(D, M, a) — Cj*| < e. 

(7.1.0) Initialize the range variables and a 2 : If Ct(D, M,ao) > tif , thena\ — (1 — fc)ao and a 2 = ao; 
otherwise, ai = ao and a 2 = (1 + /c)ao- 

(7.1.1) Use linear interpolation to predict a: Let 


a 


ai + 


c, ci (U? M, a i ) . . 

c l (D,M,a 2 )-c l (D,M,a 1 y 2 1} 


li 



(7.1.2) Check the termination condition: If 


I ci(D,M,a) -Ci | < ecf , 


then output a and terminate the algorithm. 

(7.1.3) Update ai and a 2 : 

• if a < a 1, then a 2 — au and ai — a; 

• if a > a-2, then ai — a 2 and 0:2 — a; 

• if ot\ < a < a-2 and ci(D, M,a) < c( , then au — a; 

• if ot\ < a < a.2 and ci(D , M, a) > c( , then 0:2 = a. 

(7.1.4) Start a new iteration: Go back to (7.1.1). 


7.3. Finding the Size of the Trust Region. It is well-known that the optimal objective function 
value of (7.1) is a non-increasing and piecewise affine function of 6. Therefore, we use the following search 
method based on a linear interpolation to find 8. 

Algorithm 7.2. Let 7 (d) be the optimal function value of (7.1) for any given 6. Choose an error 
tolerance e > 0 and an initial adjustment rate 0 < k < 1 . If the zero vector is not an optimal solution of 

(7.1) for 5 — d fJ with 0 < Jo < 1, then for any t > 0, the following procedure will find 0 < 5 < 1 such that 
| 7 (*) — (1 — »?)| <e. 

(7.2.0) Initialize the range variables 8\ and d 2 : If 7 (Aj) > 1 — i), then <7 = do an( l ^2 = (1 + «)do/ 
otherwise, Si — (1 — k)S 0 and d 2 = do- 

(7.2.1) Use linear interpolation to predict d: Let 


d = di + 


(1 " n) - 7(di) 
7(^2) -7(^1) 


(da -di). 


(7.2.2) Check the termination condition: If 


|7(d)- (l-»?) | < e, 


then output 8 and terminate the algorithm. 

(7.2.3) Update di and d 2 : 

• if 8 < 81, then d 2 = di and di = d; 

• if 8 > 82, then di = d 2 and d 2 = d; 

• if 8\ < 8 < 82 and 7(d) > (1 — rf), then di = d; 

• if 81 < 8 < 82 and 7 (d) < (1 — rf), then d 2 = d. 

(7.2.4) Start a new iteration: Go back to (7.2.1). 

Remark 3. In general, we use the value of d in the previous iteration of the profile optimization 
method as do- With this choice of d, k — 0.02, and e — 0.01 • i], we usually solve about 4 linear programming 
subproblems to find d in our numerical simulation runs. ■ 

7.4. Least Norm Solution of LP Subproblem. By using slack variables 1 tj and u», we can convert 
(6.1) into a linear programming problem with equality constraints and simple bound constraints: 


(7.1) min 7 subject to 

AD,A ai ,y,Ui,Vi 

0 < Ui < Ci , 0 < Vi < 1, for 1 < i < r, 

-a^k < A on < cimax - a it k for 1 < i < r, 


12 



— cri& < ADi < cf iS for 1 < i < m, 


c l (D k ,a i , k ,M i ) + {^(D k 


i M £) , A D 


dci 

+ Q^(D k ,a itk ,M i )Aa i - iq 


for 1 < i < r, 


Cd(D k ,a it k,Mi) + 



(D k ,a itk ,M i ),AD 


+ 


dcd 

da 


{D , a^ k , Af^)Acq + 


-2 for 1 < i < r. 
Pi 


The least norm solution of (7.1) can be found by solving a strictly convex quadratic programming 
problem [14] that is derived by replacing the objective function in (7.1) by 


( r m 

7 2 + J2 [«? + v i + ( Aq *) 2 ] + J2 ( A a f 

i=l i=l 

where s (— 10 -8 in our implementation) is a small positive scalar. 

Note that minimizing the objective function in (7.2) will force all v^s to be more or less the same. 
As a consequence, it produces a descent direction that gives an almost identical relative reduction of the 
drag for all the design points Mi, ... , M r . We believe that this balanced reduction of drag will prevent 
point-optimization at the sampled design points. Furthermore, the least norm solution generates a new 
design vector that stays as close to the original design vector as possible while achieving the predetermined 
percentage reduction of the drag. This is desirable in a practical design process, since if a new design deviates 
too far away from the original airfoil, then designers might have less confidence in the predicted performance 
of the new airfoil. 


8. Numerical Simulation Results. We use the NACA-0012 airfoil as the initial point for all opti- 
mization methods discussed in this section. A periodic spline representation of the NACA-0012 with 23 
control points is used to get an initial parametric model of the airfoil: 

22 22 

x = £aiPi(t)> y ~ ’^2,b i q i (t), for 0 < t < 1, 

i=0 i=0 

where pi and <"/, are spline basis functions. The shape of the airfoil can be modified by changing the values 
of the bi s. The tip and the tail of the airfoil are fixed during the optimization (i.e.. bo, &n, and 622 
are fixed) so that the chord length remains constant. As a consequence, we have 20 free-design variables 
(61 , . . . , 6 10 , 612 , . . . , 621) in the airfoil shape optimization. 

The CFD code FUN2D is used to compute the function values and gradients of the lift and drag [15]. 
Some technical details of the flow solver (for function evaluations) and the adjoint equation solver (for 
gradient evaluations) in FUN2D can be found in [1] and [16], respectively. 

We use Euler flow to demonstrate the usefulness of the profile optimization method as a robust airfoil 
optimization tool. The main reason to choose Euler flow simulation analysis is that we can complete our 
preliminary evaluation of the profile optimization method in a reasonable amount of time. 

Elliott and Peraire [6, Section 1] stated that almost any drag minimization exercise based on the Euler 
equations and applied to modern supercritical wings in cruise condition is doomed to failure. However, they 
also concluded that Euler-based optimization can be useful as a first step in the design process [7, Abstract]. 

For a given initial grid, FUN2D generates numerical approximations to the lift/drag coefficients and 
their derivatives. We use 10 -9 as the error tolerance for the flow solver in FUN2D to find an approximate 
solution of the nonlinear system derived from discretization of the Euler equation. The number of time steps 
used in solving the discrete adjoint equation is 10000 and the spring analogy method in FUN2D is used for 
grid movement [16]. 
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The initial grid for FUN2D used in this paper has 124 points around the airfoil and 32 points at the far 
field (at a distance of 50 chord lengths). The grid has 3060 nodes, 9025 faces, and 6120 elements. Fig. 8.1 
shows the initial grid around the NACA-0012. See [11] for a comparison of the numerical approximations of 
the lift and drag using this grid and other grids. 



Fig. 8.1. The grid, used for solving the Euler equation by FUN2D 

All simulation results are obtained by parallel processing of function and gradient evaluations on Coral 
(cf. Subsection 7.1). For the grid we use, it takes about 120 seconds to solve a flow problem and about 180 
seconds to solve an adjoint equation. Therefore, with 2 r processors, each iteration of the profile optimization 
method takes about 5 minutes, no matter what the number of design points is. 

We compute the least norm solution of (7.1) by the quadratic programming solver QLD, which was first 
developed by Powell [18] and then modified by Schittkowski [19]. 

Our numerical simulations are designed to examine the impact of the number of design points and 
the target lift value on the profile optimization method. Also, we shall compare the simulation results 
obtained by the profile optimization method with the multipoint optimization method and the expected 
value optimization method. Therefore, we consider the following cases for numerical simulation: 

1. the profile optimization with 4 design points equally spaced in the Mach range [0.7, 0.8], c\ — 0.4, a 
fixed percentage reduction rate rj — 5%, and a termination parameter e = 0.3%; 

2. the profile optimization with 3 design points equally spaced in the Mach range [0.7, 0.8], c ; * = 0.4, a 
fixed percentage reduction rate rj — 5%, and a termination parameter e = 0.3%; 

3. the profile optimization with 8 design points equally spaced in the Mach range [0.7, 0.8], c/ = 0.4, a 
fixed percentage reduction rate r/ — 5%, and a termination parameter e — 0.3%; 

4. the profile optimization with 4 design points equally spaced in the Mach range [0.7, 0.8], c ; * = 0.2, a 
fixed percentage reduction rate rj — 5%, and a termination parameter e = 0.3%; 

5. the multipoint airfoil optimization with 4 design points equally spaced in the Mach range [0.7, 0.8], 

W\ — = 1/6, W 2 — Wz — 1/3, and c ( * — 0.4; 

6. the expected value airfoil optimization with adaptive changes of w* corresponding to 4 randomized 
integration points, and c ( * = 0.4. 

The first case is the benchmark case. The next two cases are included to examine the impact of the 
number of design points on the profile optimization. We use Case 4 to examine the impact of the target 
lift value on the profile optimization. The last two cases are for comparison of the profile optimization with 
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other optimization methods. The profiles of the drag coefficient for a given c ; * are plotted by using 24 equally 
spaced Mach numbers in [0.7, 0.8]. 

8.1. Impact of the Number of Design Points. The profile optimization method terminates after 
69, 67, and 57 iterations for 3, 4, and 8 design points, respectively. The number of design points has no 
significant impact on the optimal airfoils generated by the profile optimization method. All optimal airfoils 
have similar shapes (cf. Table 8.1) and only a marginal drag reduction is achieved by adding more iterations 
(cf. Table 8.2). 


Table 8.1 

The relative differences ( in percentage ) of the design vectors 



NACA-0012 

3 Points 

4 Points 

8 Points 

NACA-0012 

- 

24.5 

25.7 

23.8 

3 Points 

24.7 

- 

1.7 

1.2 

4 Points 

25.9 

1.8 

- 

2.6 

8 Points 

23.9 

1.2 

2.6 

- 


Each number in Table 8.1 denotes the relative difference of two design vectors: ||D r — D c ||/||D r ||, where 
D r (or D c ) is either the NACA-0012 or the optimal design vector obtained by the profile optimization 
method with the number of design points listed in the corresponding row (or column). 

A typical optimal airfoil is given in Fig. 8.2(b). This Mach wave plot illustrates how one strong shock 
wave for the NACA-0012 is reduced to several weaker shock waves for the optimal airfoil. 



(a) NACA-0012 


(b) Profile Optimization 


Fig. 8.2. The Mach waves (for the Mach number 0.796) around the NACA-0012 and an optimal airfoil generated by the 
profile optimization method 


With 4 design points, the profile optimization method achieves consistent drag reduction over the Mach 
range [0.7, 0.8] (cf. Fig. 8.3(a) and 8.4(a)). There is no random distortion of airfoil shapes during the 
optimization process (cf. Fig. 8.3(b) and 8.4(b)). 

It is worth pointing out that the iterates generated by the profile optimization method are also indepen- 
dent of the number of design points. For example, the 56-th iterate generated by the profile optimization 
method for 4 design points differs from the 56-th iterate corresponding to 3 (or 8) design points by 0.35% 
(or 1.2%). The profiles of the drag coefficient for optimal airfoils generated by the profile optimization with 
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3 and 4 design points, respectively, are almost identical. The profile of the drag coefficient for the optimal 
airfoil corresponding to 8 design points is very similar to the drag profile of the 57-th iterate generated by 
the profile optimization with 4 design points. 

8.2. Impact of the Target Lift Value. The target lift value c ( * has a quite significant impact on the 
optimal airfoil generated by the profile optimization method. 


14R 
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(a) Profile Optimization With Lift at 0.2 


Design Points 
NACA-001 2 
the 6-th iterate 
the 1 2-th iterate 
the 1 8-th iterate 
the 24-th iterate 
the 33-th iterate 



0.7 0.71 0.72 0.73 


0.74 0.75 0.76 0.77 
Mach Number 


0.78 0.79 0.8 


(b) Profile Optimization With Lift at 0.2 



Fig. 8.3. The two figures illustrate the behaviors of the iterates generated by the profile optimization method, with cj = 0.2 
and 4 equally spaced design points. Fig. 8.3(a) shows profiles of the ding coefficient versus the Mach number for the iterates 
and Fig. 8.3(b) shows the shapes of airfoils corresponding to the iterates. (Note. The airfoils are shifted to improve the 
visibility.) 


With Ci — 0.2, the NACA-0012 is almost optimal for low Mach numbers. Therefore, we can not 
simultaneously reduce the drag over the Mach range [0.7, 0.8]. The profile optimization method tries to keep 
the drag as low as possible for Mach numbers between 0.7 and 0.75, while reducing the drag significantly for 
Mach numbers between 0.75 and 0.8 (cf. Fig. 8.3(a)). 

Note that the drag bucket («.e., a drop of the drag before its dramatic rise) occurs near M — 0.73 for 
the NACA-0012 and the drag bucket for the optimal airfoil occurs near M — 0.78. Such a delay of the drag 
rise is desirable. The only compromise made by the optimal airfoil to the NACA-0012 is a small increase of 
the drag around the original drag bucket for the NACA-0012, which is quite reasonable. 

For the higher target lift c ( * = 0.4, the optimal airfoil shape (shown in Fig. 8.4(b)) deviates more from 
the NACA-0012. Also, the 33rd iterate generated by the profile optimization method for c ; * = 0.4 differs 
from the 33rd iterate for c ; * = 0.2 by 10%. 

Note that the shape variations of the airfoils corresponding to iterates generated by the profile optimiza- 
tion method follow a similar trend no matter what the target lift is: the aft end thickens while the front 
section narrows. This pushes the shock location towards the aft region of the airfoil. 

A similar airfoil with thin front and fat tail was obtained by Elliott and Peraire [6, Figure 34] in a 
totally different context, where they used two thickness design variables and one camber design variable for 
a lift-constrained drag optimization with 2-D separated viscous flow and an additional area constraint. The 
initial airfoil used by them is also the NACA-0012. Their explanation of the merit of an airfoil with thinner 
front and fatter tail over the NACA-0012 is that the thickness distribution has been redistributed such that 
the maximum is further aft (i.e., closer to the tail), like the early natural laminar flow airfoils, this delays 
the start of the adverse pressure gradient to aft of that maximum thickness point. 
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8.3. Impact of the Optimization Strategy. From the analysis given in Sections 3 and 5, we know 
that if the design points and weights are fixed as in either (5.1) or (5.3), then a severe degradation in the 
off-design performance is likely when the number of design points (i.e.. 4) is much smaller than the number 
of free-design variables (i.e... 20). Fig. 8.5(a) clearly reveals the point-optimization features of the optimal 
airfoil generated by the multipoint optimization method. Fig. 8.4(c) suggests that the point-optimization 
behavior does not occur until the end of the optimization iterations, since the drag curve for the 25-th iterate 
does not have any drag trough. 

The expected value optimization method [12] adds a random perturbation to the integration points 
from iteration to iteration and adjusts the weights accordingly so that a severe degradation in the off-design 
performance can be avoided. Fig. 8.4(e) shows that this technique allows the optimizer to reduce the drag 
while avoiding point-optimization. 

The profile optimization avoids the problem of point-optimization at the design points (cf. Fig. 8.3(a) 
and 8.4(a)) by using a smart descent direction to achieve consistent drag reduction over the whole Mach 
range. 

Table 8.2 

Relative drag reduction rates (in percentage) 



Max 

Min 

Average 

Profile Optimization (3 Points) 

82.9 

5.1 

70.2 

Profile Optimization (4 Points) 

82.7 

1.4 

69.7 

Profile Optimization (8 Points) 

82.6 

4.1 

67.9 

Multipoint Optimization (4 Points) 

89.4 

46.2 

81.1 

Expected Value Optimization (4 Points) 

89.4 

49.6 

82.6 


Table 8.2 gives an overall sense of how each optimization method reduces the drag coefficient with 
c*i — 0.4. Let Mi,M 2 ,..., M 2 4 be 24 equally spaced points in [0.7, 0.8]. For each M iy let Cd,i (or c^,*) be 
the drag coefficient of the NACA-0012 (or an optimal airfoil) with the lift coefficient fixed at 0.4. Then the 
relative reduction of the drag at each Mach number M t is (cd,i — Ci^)/c d l . The numbers in “Max”, “Min”, 
and “Average” columns of Table 8.2 correspond to 


max 

1<<<24 


c d,i Cd,i 

Cd,i 


c d,i &d,i 1 

min — ! -, and 

1 


/ 24 24 \ / / 24 \ 

( ^2 c d,i - ^2 Cd,i ) / I ^2 Cd,i I , respectively. 

Note also that the expected value optimization method (as well as the multipoint optimization method) 
resulted in an optimal airfoil with a drag curve that is consistently lower over the Mach range than the 
profile optimization method (cf. Fig. 8.5(a)). It seems that it should be possible for the profile optimization 
method to find a descent direction for a consistent drag reduction over the whole Mach range [0.7, 0.8], by 
moving toward the optimal solution given by either the multipoint optimization method or the expected 
value optimization method. However, due to nonlinearity of the drag coefficient (with respect to the design 
variables), the profile optimization method can not find such a descent direction by using the local derivative 
information only. 


9. Conclusion. In this paper, we introduce a new robust optimization scheme called the profile opti- 
mization method and use airfoil optimization under uncertain flight conditions as a case study to evaluate 
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(a) Profile Optimization 


(b) Profile Optimization 


0.025 


0.02 


0.015 


0.01 


0.005 


0 

Design Points 


NACA-001 2 

* 

the 13-th iterate 

□ 

the 26-th iterate 

> 

the 39-th iterate 

0 

the 52-th iterate 

<1 

the 69-th iterate 



0.7 0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.8 
(c) Multipoint Optimization 



(e) Expected Value Optimization 




From Left to Right: NACA-001 2, 1 3-th, 26-th, 39-th, 52-th, and 69-th Iterates 
(d) Multipoint Optimization 
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Fig. 8.4. These figures show the behaviors of iterates generated by various optimization methods with cj = 0.4 and 4 
design points. 


this method. We used an Euler-based CFD code, which does not include viscous effects, to test the profile 
optimization method. Because of this lack of fidelity, the generated airfoils may be somewhat unrealistic. 

The profile optimization method adaptively adjusts the weights in a minimax optimization formulation 
to find a drag reduction direction for all design conditions, which leads to a consistent drag reduction over 
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(b) Optimal Airfoils 



From Left to Right: Profile Opt, Multipoint Opt, and Expected Value Opt 


Fig. 8.5. Fig. 8.5(a) compares the profiles of the drag coefficient for various optimal airfoils and Fig. 8.5(b) compares 
the optimal airfoil shapes. 


a given range of Mach numbers from iteration to iteration. 

Without adaptive adjustment of weights, we prove that it is necessary to use at least (m + 1) design 
points, where m is the number of free-design variables, in order to avoid point-optimization at selected design 
points. 

Our numerical results demonstrate that the profile optimization method is not sensitive to the number of 
selected design points. With 20 free geometric design variables and as little as 4 design points, the optimal 
airfoil generated by the profile optimization method is smooth and has no degradation in the off-design 
performance. 

The profile optimization method can be easily modified to solve other optimization problems under 
uncertainty by replacing c<j with another performance measurement function and M with other uncertain 
parameters. 

The profile optimization method also has the potential of becoming a practical design tool for optimiza- 
tion under uncertainty. The use of a small number of sampled design points from the range of uncertain 
parameters makes the profile optimization method computationally affordable. The consistent performance 
improvements over the range of uncertain parameters from iteration to iteration allows a designer to stop 
the iterative process at any time with an improved new design. 
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